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ABSTRACT 

In This paper we consider elliptic and hyperbolic problems in unbounded 
regions. These problems, when one wants to solve them numerically, have the 
difficulty of prescribing boundary conditions at infinity. Computationally one needs a 
finite region in which to solve these problems. The corresponding conditions at 
infinity imposed on the finite distance boundaries should dictate the boundary condi- 
tions at infinity and be accurate with respect to the interior numerical scheme. Such 
boundary conditions are commonly referred to as absorbing boundary conditions. This 
paper presents a survey and cover our own treatment on these boundary conditions 
for wave-like equations 
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ABSORBING BOUNDARY CONDITIONS 
FOR EXTERIOR PROBLEMS 

by 

S. [. Hariharan* 


l • Intr oduct ion 

Many formulations arising from physical nature yield problems in un- 
bounded regions. Mathematical formulations of such problems yield govern- 
ing partial differential equations in or near a given domain in such a fashion 
that: i) the equations may be linear but with non-constant coefficients, or 
ii) the equations may be nonlinear, but at large distances essentially be- 
have linearly and with constant coefficients. This note presents a survey 
of the treatment of such problems, when the desired solutions are governed 
by elliptic or hyperbolic partial differential equations. These problems are 
called exterior problems and commonly arise in the fields of aerodynamics, 
meteorology, electromagnetic scattering, and atmospheric acoustical wave 
propagation. The main difficulties with these problems are the boundary 
conditions that need to be prescribed at large distances from the region of 
interest. Usually only an asymptotic behavior is known. Such conditions 
may be sufficient to check the well- posed ness of the problem, however, if 
one wants to compute the solutions of these problems numerically, infinite 
distances need to be truncated to finite distances. The boundary conditions 
imposed on these finite distance boundaries should dictate the behavior at 
infinity and be accurate with the interior numerical scheme. Furthermore, 
the shorter the distances, the more efficient the solutions in terms of com- 
putation time required. This consideration is the essential need in several 
problems, depending on the problem and the kind of computer that is used. 
Typical of such cases are most three dimensional problems. The intention 
here is to present some available techniques to overcome this difficulty, in- 
cluding application to some model problems. 

To begin, let us illustrate some concepts using simple one dimensional 
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model problems. First, let us consider an elliptic problem where the field 
equations are reduced wave equations. 

Consider a slab of thickness L with a varying index of refraction n(x). 
Left of this slab, let the media be homogeneous with index of refraction l 
and right of the slab (x > L) let the index of refraction be n 0 (> l), which 
is a constant. Such problems are common in optics as well as geophysical 
waves. This situation is illustrated in Figure 6.1 The governing equations 
are as follows: 



Figure 1 One dimensional model 


Then (1.1) needs to be solved in [0, L '\ , [U > L) say with boundary condi- 
tions 


a,u(0) = g 

(1-2) 

B 2 a(L') =0. 


a, a 1 continents on interfaces 

(l.i) 


where D\ and B 2 are boundary operators and g are given are to be chosen 
according to the physics of the problem. For example, if there is a unit 
*> 



amplitude wave traveling from — oo incident on the slab, then for x < 0 the 
solution may be written as 

«*(x) = e lkx + R{k)e - tkx (1.5) 

where R(k) is the reflect ion coefficient and measures the part of the wave 
which is reflected from the boundary x = 0. Eliminating R(k ) from (1.5) 
yields 

u'(0) + iku( 0) = 2 ik (1.6) 

which is called an inflow condition. Thus in view of (1.2) B t = d/dx 4- 
ik and g = 2 ik. Now we can do a similar calculation to obtain a boundary 
condition at L'. For x > L all the waves transmit and do not reflect back as 
there are no other boundaries for x > L. In this region the solution may be 
written as 

u(x) = T(k)e tknnx (1.7) 

where T(k) is the transmission coefficient which measures the transmitted 
part of the wave. Eliminating T(k) in (1.7) we obtain 

n'{L') — iknou(L') = 0 (1.8) 

Again comparing with (1.3) we see that the operator B 2 is 

#2 = t — ikn 0 (1.9) 

ax 

The boundary condition (1.8) is the desired absorbing boundary condition 
for this problem, which is exact. 

We note that this condition could have been applied exactly at x = L. 
and this, as we will see. does not always hold in higher dimensions. There is 
another technical difficulty in this type of problem. One can find that even 
if n 2 (x) = constant there is a countable set of wave numbers {k Tl \ called 
the interior resonant values for which the continuous problem does not have 
solutions. When one deals with higher dimensional problems it is not easy 
to calculate such frequencies when one has a general shaped region where 
the solution is sought. 
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Figure 2 One dimensional model 


As an example of hyperbolic problems and associated boundary condi- 
tions, let us consider again a one dimensional model. In particular let us 
consider the following problem (see figure 2): 


c 2 {x) 


it — ^ix 


2 ^tt — 
c 0 


0 < x < L 


x > L 


u, u x , u t continuous on x = L 
tt(0,i) = f(t) 

Hu(L',l) = 0 (£/>£,). 


( 1 . 10 ) 

( 111 ) 
( 1 . 12 ) 
(1 .13) 


This problem tells us that the wave is propagating from a source at i = 0 
and the wave speed changes in 0 < x < L and for x > L the waves do not 
reliect in this region. Equation (1.13) dictates a no reflection condition in 
this region and D is an operator which is to be determined to fit this physical 
nature. For x > L the outgoing wave solution can be written as 


Which yields 


u[x. t) = — c 0 t) 


c'o + ut — 0 


Thus the absorbing condition to be imposed at x — L' is (1.15) and 


„ _ 0 d 

11 = c °ox + al 


(in) 


(1.13) 


(1.16) 
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We see from both one dimensional examples that the absorbing bound- 
ary conditions are exact and easy to obtain. The difficulties arise in higher 
dimensions as we will see. Nevertheless the first problem we described here 
serves as a model for inverse problems that can be investigated numerically 
to calculate n 2 (x), if the no reflection condition B 2 U = 0 is given. This is 
described in Dunn and llariharanl 6 L 

The approach taken is as follows: fn section two we describe a local 
boundary condition procedure for exterior elliptic problems. In section three 
we describe a nonlocal boundary condition for similar problems. In section 
four we describe a new method that combines the previous two procedures 
but more accurate than either. Finally, in section five we describe the meth- 
ods known for time dependent problems and describe their use for a practical 
problem. 

2 . Loc a I B ound ar y Conditions Procedure 

We pointed out in the last section that the absorbing boundary condi- 
tions are exact in one dimension and do not hold in higher dimensions. We 
want to investigate this statement a little further and provide a treatment for 
this, that of Bayliss, Grunzburger and TurkelJ 2 - For illustrative purposes let 
us consider the following problem that commonly arises in two dimensional 
acoustic scattering. 

Let H be a simply connected, bounded domain with boundary T and its 
exterior H Then the problem is (see figure 3): 

Ait -f k 2 u = 0 in Q (2.1) 

^ = g on r (2.2) 

an 

u satisfies a radiation condition (2.3) 

Briefly this problem has a physical meaning that, there is a wave incident on 
n whose boundary T is a perfect reflector and the reflected waves all decay at 
infinity, u measures the scattered part of the wave and g is the contribution 
from the incident wave. Regardless of the boundary condition on T, the issue 
here is the radiation condition. Note that the radiation condition is another 



Figure 3 Two dimensional scattering problem 


name for the absorbing condition. It is known that at large distances the 
solution of (2.L) behaves like 


+ h(fl + M?J + .. 

r r 2 

( 2 . 1 ) 

The first part of this solution dictates outgoing waves and the second part 
dictates incoming waves, which can easily be seen from the signs e Tlkr . Thus 
in order for the radiation condition to be satisfied we must pick only the first 
part. 


u ~ 


■i kr 


vA 


a o {0) + 


ai(0) a 3 (0) 


+ 


+ 


,-ikr 


V r 


bo(O) 


i.e. 


a 



^o(^) + 


a y { 9 ) 0,2(9) 


+ 



(2.5) 


Comparing this to the one dimensional case in (1.7), shows that instead of 

a single transmission coefficient we have many coefficients a t ,i - 0, 1.2 

Thus the one dimensional way of eliminating the transmission coefficient 
does not exactly work in this situation. This shows the difficulty in higher 
dimensions. However, if we settle for loss of accuracy, in particular eliminat- 
ing clq(9) we have 


>ikr 


u r — iku = — ( (Iq ( $ ) + 

2r 2 r 

e lkr ( ay( 9 ) ■ 2 a 2 ( 9 ) 


«i(0) , a 2 (0) 


+ + • • •) 


+ 


+ 


Thus we see that 
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s/r \ r* r° 

r'i(u r — iku) = 0 ( l/r) 


( 2 . 6 ) 


C2.7) 


or simply n (u r — iku) = 0 which is Sommerfeld’s radiation condition 

in two dimensions. 

Taking a closer look at the right hand side of (2.6), in particular, the 
first expression, we see that it is nothing more than —u/2r. Thus we have 


u r — iku -f ^-u — 0 ^1/r 5 / 2 ^ (2-8) 

Comparing (2.7) and (2.8) shows the error introduced in the boundary con- 
dition at large distances is reduced by a factor of 1/r. This was first studied 
in connection with numerical implementation for time dependent problems 
by Bayliss and Turkell 3 ! and also by Engquist and Majdaf 7 !. For time har- 
monic cases with time dependence of the form e~ lkt , both of their boundary 
conditions reduce to the form (2.8). Bayliss, Gunzburger and Turkel^ gen- 
eralized these conditions. They defined the boundary operator on the left 
side of (2.8) by 

B { = — -ik+ l/2r (2.9) 

dr 

When (2.8) is explicitly written it has the form 


,ikr 


B\u = 




ax(9) 2oj(9) 

o T o 


+ 


( 2 . 10 ) 


They observed that now the second coefficient ai(0) can also be eliminated. 
If one does that by setting v = Byu it is readily verified that 

dv 


37 - ikv + - r v = 0 ('A 9 ' 2 ) 


or 


= {i-‘ k+ Tr) [Fr~ ik+ Tr) “ = ° ( l/r< " 2 ) 


d 


( 2 . 11 ) 


This process can be repeated to generalize 


where 


i=i x 

(2.13) 

Thus in order to implement this approach the problem can 
a truncated region (figure (4)) as follows: 

be considered in 

Seek w such that 


Aw + k 2 w = 0 in Or 

(2.11) 

dw 

a7 =5 ° nl 

(2.15) 

B rn w = 0 on Toq. 

(2.16) 



Figure 4 Computational domain 


Where is a circle where the approximate boundary conditions are 
to be imposed. 

In the region Or one can use a finite element formulation of the problem 
similar to the one used in the next section to seek the solution. But there 
are some difficulties. 


First, the governing equation is a second order differential equation. 
Even the second order boundary condition involves second order derivatives 
in the radial direction. And all the higher order conditions (2. 13) have higher 
order radial operators. This can be overcome by appealing to the differential 
equation itself. Rewriting the equation in cylindrical coordinates 


(A + k 2 )w 


d 2 w 1 dw 1 d 2 w 

< )r 2 r dr r 2 36 2 


+ k 2 w — 0 


(2.17) 
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From (2.17) we see that second order radial operators can be translated into 
single First order radial derivatives plus derivatives in the tangential direction. 

The second difficulty is the following. In comparing (2.1) - (2.3) and 
(2.14) - (2.16) we see that there is an a priori error introduced (irrespective 
of the numerical scheme used) of the order 0(1 //? im+l / 2 ), where R is the 
distance of F^ from the origin. This of course depends on the order of the 
operator too. The question is what can we say about || u — w || in a suitable 
norm. The answer is not available in two dimensions. However, Bayliss, 
Gunzburger and Turkel were able to prove for corresponding operators B rn 
in three dimensions for m = l and m = 2 the following theorem. Let 
w = u — w, then 


ll®ll(r )=C/r m+l , (2.18) 

where C depends on k and Too and the norm is defined on the surface as 

II <0 11(0= l I «> I 2 dA 

This theorem tells that the error on the artificial boundary in that given norm 
is inversely proportional to l/r m+l and this error remains the same in the 
interior boundary too. The key point here is that this bound is dominated by 
C which depends on the wave number k. For three dimensions as k — > 0 this 
estimate is still valid. However, in two dimensions, even though it is possible 
to get such an estimate as the wave number becomes smaller, this constant 
can grow larger. In fact, ill two dimensions there is a logarithmic branch 
point as k — 0 and the constant C can grow very large. For this reason, for 
low frequency cases the problem must be examined very carefully. Such a 
treatment is technical and it is available in Hariharan and MacCamy l2 i. 

These difficulties are pertinent to the reduced wave equations. However, 
other strongly elliptic cases can be handled without difficulties. For example 
if we consider Laplace's equation Au = 0, then the solution at infinity either 
has to be bounded or behaves logarithmically with a given behavior. To 
apply the above process let us consider the following problem: 

Au = 0 in fl + 


(2,19) 
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a = (j on r 

u, ~ log r as r = \/a; 2 + y 2 


( 2 . 20 ) 

( 2 . 21 ) 


— ► oo. 


Solutions for this problem can be written as 

OO 

u(r,0) — logr+ ^ ^ cosnO (2.22) 

H~0 

This can be used to obtain boundary conditions on the artificial boundary 
r^. Set \) = u — logr — a 0 . Eliminating a t we see that. 


dv v . i . 


Similarly eliminating «2 we obtain 


(2.2.3) 



( 2 . 2 . 1 ) 


This process can be repeated to obtain higher order conditions. The constant 
a 0 can be obtained by averaging the solution on Poo. In section 4 we provide 
an alternative treatment to resolve the difficulties. In fact we implemented 
the conditions (2.23) and (2.21) along with ours to compare the efficiency of 
our procedure. 


3. Nonlocal Bou ndary Condi tions Procedure 


This procedure is a little more difficult to handle than the previous 
one. Hut, it does not introduce an a priori error due to placing the arti- 
ficial boundary at finite distances. There are two different versions of this 
procedure. The first one originated in the work of Fix and Marin— which 
was done for situations of wave guides and was made general in two dimen- 
sions by MacCamy and Marini 15 !. The second version is due to .Johnson and 
Nedelec 1 tj ! , which was done independently but has similarities in the ap- 
proach. Extension to three dimensions was done by Aziz and Kellogg 1 '. All 
these procedures are done in view of implementing finite element methods. 
As a result, a reasonable amount of analysis is available for this method. It is 
almost impossible to summarize all of it here. We choose the method of Mac- 
Camy and Marin and describe how the nonlocal conditions are treated. The 


to 



Figure 6 Two dimensional interface problem 


attractive feature of this method is that it is appropriate to interface prob- 
lems analogous to the one dimensional model that we described in section 
one. Let us present the model problem in two dimensions. Again this is a 
scattering type of a problem, but this time the scatterer is an inhomogeneous 
penetrable body. The problem is as follows: 

A u + k 2 ri 2 (x)u = / in Q (3-1) 

Aw 4- k 2 u = 0 in il + (3.2) 


a 


= u + on r 


du 

dn 


du + 

dn 


on r 


(3.3) 

(3.4) 


u — iii satisfies Sommerfeld’s radiation at infinity. (3.5) 

equations (3.1) and (3.2) say that the index of refraction of the media in Q is 
n 2 (x),x = [x,y) and outside is just a constant. In (3.5) u,(x) is a prescribed 
incident field which satisfies (3.2). (3.3) and (3.1) are continuity conditions 
of the solution on l\ In many applications they may not be continuous, but 
rather may have jumps. These jumps can be treated without much difficulty. 
To describe the procedure let us extend n 2 (x) in Z 2 so that 



in n 
in n + . 


and extend / in a similar manner, so that 


in n 
in + 


II 



So that the governing equations will be 

Au + k 2 nu = g 


in Z 2 . 


(3.6) 


To do numerical calculations let us truncate the infinite region by a circle 
Too into a finite one as depicted in Figure 6. Denote the truncated region 
by fir- 



A, 


Figure 6 Computational domain 


A standard way of handling this problem is to use variational methods. 
Suppose u is a solution and v is an arbitrary differentiable function in fir- 
If we multiply (3.6) by v , integrate over fir and use integration by parts we 
find 

Va-Vu + fc 2 / n uv + / v— = / Q v (3-7) 

T ./ n T ” n Jn T 

This variational form is the basis of a numerical method. But, in order 
to implement this numerically, one should have enough information about 
Ou/dn on r^. This in fact requires the knowledge of the solution exterior 
to Too in the region Aoq. Recalling that in this region u satisfies 


-I 


Au + k 2 u = 0 


This suggests that one rnav seek, an integral representation of the solu- 
tions. In particular let us consider 


where Gk is the free space Green’s function and has the form 


G k {x,y) = - - H^ x) {k |x-y|) (3.9) 

The representation (3.8) is known as the simple layer representation, which 
satisfies the Helmholtz equation and satisfies the radiation condition at in- 
finity. These two facts are easily verified by applying the Helmholtz operator 
to (3.8) and by realizing the asymptotic behavior of (3.9) for large values of 
| x | which yields the form (2.5). In order to calculate the normal deriva- 
tive of u on Too, (3.8) together with the standard jump relation in potential 
theory yields 


~(x)=^(x)+/ "(y)J ; G t (x,y)o!» s 


du t (x) „ 

+ x € r„ ( 3 . 10) 


Thus ()a/()n can be calculated once er is known. One can use (3.8) to obtain 
a singular integral equation of the first kind to determine a when x is on 


r«. 


u s = u - Ui= / cr(y)Gfc(x,y)ds y , x G T c 

Jr™ 

Let us denote this equation in an operator form 


(3. 11) 


u„(x) = Gk o cr(x), x € Too (3.12) 

Then a can be formally inverted to obtain its value of the form 

cr = G^ l ou a (3.13) 

But there is a technical difficulty here. G^ 1 exists provided - k 2 is not an 
eigenvalue of the operator A with zero Dirichlet conditi(?n on Too. This result 
is familiar in the diffraction theory and it is analogous to the statement that 
we made for the corresponding one dimensional problem. That is to say, 
such values of k will be the interior resonant values of Qt- But there is 
a way to treat this problem. Ursell^ l8 l proposed to seek solutions in the 
exterior not only by a simple ltiyer operator as in (3.8), but combined with 
a double layer representation. This slightly complicates the explanation of 
the present method, but nevertheless, can be done. For the moment we shall 
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assume k is not a resonant value so that (3.13) holds. Then (3.10) takes the 
form (for x 6 Too): 

^-1 , v dG k du,(x) n 

Gfc OMs ^ a^T^ x,y ^ 5y + ~d^~ ^ 3 ‘ 14 ^ 

The right hand side of (3.M) is another integral operator acting on a, so 
that (3.11) can be written (for x £ Too) as: 

i)uf <)n = T k (u) - T k (u,) - (3. 15) 

Thus, the normal derivative of u is a functional of u on Too. A closer look 
at (3.M), shows that to calculate du/dn at each point on r^, we need the 
knowledge of a on the entire boundary Too. Thus the boundary condition is 
nonlocal. Properties of the operator T k are discussed in detail in reference! 15 !. 
In particular it is a bounded linear operator. 

Returning to the variational formulation, (3.7) together with (3.15) one 

has 

a(ii.v) = — / Vu-Vv + k 2 fitiv T / >' T k {u) = 

J O r Q t *M' jo 

I 9 V + [ {Tk(ui) + -r-^)v = t (*-*) (3.16) 

Jn T Jr c* dn 

To provide a brief numerical implementation of (3.16) we seek approximate 
solutions u h such that 

a(u h k v h ) = F(v h ) (3.17) 

for all i’h £ S h . 

where u h and n/ t are in a finite dimensional subspace S h of an infinite di- 
mensional space $ where a is sought. Then (3.17) is made equivalent to a 
matrix problem by selecting a basis {pi, p 2 < • • • • py } for S h . This says a 
can be approximated by 

y 

= Y2 q y p i ( 3 - 18 ) 

j - » 

which satisfies 



It 




(3.19) 


Then (3.19) is the matrix problem 


k q = f 


(3.20) 


where q = ( 71,721 , 7 /v) T is the vector of weights in (3.20) and f is the 

source term. 

k L = = - / Vpy • v<p t + k 2 / nVjPi + / <piTk(<Pj) 

■'Dr ^ n r '/roc 

(3.21) 

We can now describe how the nonlocal boundary condition (3.15) is used. 
We see from (3.21) that calculating (3.15) is equivalent to computing the 
integrals 


I <p % Tk(<Pj)ds 
Jr 00 


(3.22) 


for the basis functions pi,p 2 , •• •<pjv of the approximation space S h . This 
computation is carried out in a straightforward manner using (3.12) and 
(3. 11). First solve for cr t , i — !, • • •• N. 


I, 


ffi(y)Gk(jt,y)* y = p,(x) , x s r s 


(3-23) 


and compute T k (p l )(x) (for x G Too) from 

T k {p t ){x) = ^(x) + ^ (r l {y)^C k {x,y)dsy (3.24) 

There are effective procedures to implement (3.23) and (3.24) in two 
dimensions which can be found in MacCamy and Marin 115 !, and also in 
conjunction with an integral equation treatment to this problem, when n 2 (x) 
is a constant, found in Hariharan and MacCamy! 121 . 

Some final remarks are needed. We assumed u and du/dn are continu- 
ous on the boundary f. If they have jump discontinuities of the form 


The variational form should be modified. These modifications can be found 
in reference! 15 ! and in Bielak, MacCamy and McGhee! 4 ! together with nu- 
merical implementation. 

4. Infinit e or der radiati on condition procedure. 

In this section we provide a description of a new method due to Canuto, 
Hariharan and Lustman! 5 !. From the discussions of local boundary conditions 
procedure, we saw that errors introduced in computations are twofold. The 
first one is the error which varies according to the farfield distance and the 
order of the boundary condition used. The second one is due to implementa- 
tion of finite element method. In the nonlocal boundary condition procedure 
we see that error is essentially only due to the finite element method and 
the solution of the integral equation on the boundary No other error is 
introduced. The local conditions can be improved by either increasing the 
distance of the artificial boundary or increasing the order of the operator 
B rn so that the dominating error will be essentially due to the finite element 
implementation. Error estimates, in such procedures are usually at the best 
of 0 (h 2 ) for this type of problem. A question which then can be asked is can 
we improve the accuracy of the solutions by a method which does not have 
an a priori error due to placing the artificial boundary and at the same time 
achieve error in computation 0 (h M ) , where M is the number of elements 
used. It is possible, if one uses spectral methods, and a brief theoretical 
discussion of obtaining such accuracy is found in Eustman! 14 ! for simple one 
dimensional problems. The method we will describe here does not have any 
theoretical error estimates yet. At this point we only have numerical evi- 
dences. Theoretical error estimates are available oniy in one dimension, and 
extension to higher dimensions is still much in need of work. The proce- 
dure we are going to describe here works for a general second order elliptic 
problem, exterior to a given domain in two dimensions. The structure is as 
follows: 


Lu = f in Q + (4.1) 

u — g on T (4.2) 

Boon = 0 as | x |= \fx? + y 2 — ► oo (4.3) 


Where H + and V are same as that described in section 2. B 0 0 is a boundary 
16 


operator at oo. For illustrative purposes and to appraise numerical results 
let us consider the problem given through (2.19) - ( 2 . 21 ), which we repeat 
here. 


A a =0 in fi + (4.1) 

u = g on T (4.5) 

u ~ log r as r = s/x^ + y 2 — ► oo (4-6) 

Again the goal here is treat the condition (4.6) numerically. Solution in the 
farfield may be sought through separation of variables of the form 

,<(r, v) = log r + Y, H- 7 ) 

k 


where (r, p) are the polar coordinates in the plane. Note that the right hand 
side of (4.8) satisfies the radiation condition (4.6) at infinity. The coefficients 
ak are unknown. The approach here consists of expressing each coefficient 
a k as a functional of u, rather than eliminating a finite number of them using 
differential operations described in section two. From (4.7) we see that for 
any r > 0 , rifc/r^ is the fc-th Fourier coefficient of the periodic function of 
<p h-+ u(r, p). Thus we can invert the coefficient to obtain 

7 n = , 7 - / u{r,d)e~' kd d6 = u k {r) (4.8) 

H fc I 2 tt J o 

If we differentiate (4.7) with respect to r 


“r{r,p) = 

r 


i-E 1*1 

k 


r \k\ c 


(4.9) 


and use (18), we obtain an integro-difFerential relation on circle of radius r, 


as follows: 



r 


| k | e lk{ ' p ~ 9) u{r,9)d9 

(4.10) 

- K o u 

(4.11) 

r 
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where K o u is the convolution of u with singular kernel 

1 °° 

K(rj) = ~ ^ fccos krj (4.12) 

7T 

fc = 1 

Again we consider a truncated region Qt with an artificial boundary 
Too where the condition corresponding to that of (4. LI) will be imposed. 
Note at this point that (4.11) is similar to the condition (3.15) except it 
is imposed on a circle of radius r = R <x . It is not a necessary restriction 
that I’oo should be a circle. Llut the restriction guarantees accuracy, if the 
approximate solution is a trigonometric polynomial on Too, which is the case 
if a spectral Fourier method is used in the angular direction. Suppose the 
approximate solution on F^ to be u' v , which is a trigonometric polynomial 
of degree .'V : 

«"(Koo,*>) = \ k I («») c' kV ( I I I) 

\k\<N 

The asterisk in the summation indicates periodicity in the <p direction (i.e., 
u ( R. 00 ) = n{Roo))- Then (4.16) says that the integral operator K pro- 
duces a new polynomial of degree A r , whose Fourier coefficients are obtained 
from those of u^IRqo, •) by multiplication by the modulus of the wave num- 
ber | k | . If ir v is known on Foo through its values u N (R <x> , <Pj) at equally- 
spaced points jn/iX , j = 0, !,••••, 2 N — 1 , then K o u N can be computed 
at the same nodes exactly and efficiently, first, by Fourier transforming the 
values of ir v to get its coefficients, then multiplying by the modulus of the 
wave numbers and finally by Fourier transforming back to get the point 
values of K o u ,v . This takes order of iVlog 2 A r operations if one uses fast 
Fourier transforms, which are always used in this type of calculations. Thus 
the spectral solution is required to satisfy the radiation condition 

a; v = -j- [1 - K o «"] on Too (4.15) 

/loo 

Unlike the fatnily of conditions described in section 2, this is precisely the 
same boundary condition satisfied by the exact solution except for the trun- 
cation error which comes from using a finite number of modes. No other error 
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So that at s = — 1 , r( — 1, 0) = R(9). Demanding r( 1, 0) = we determine 
the stretching parameter a(0) as 

<*(«)=£ l°g(/W«(0)) (4.19) 

This transformation in turn changes the Laplace equation into the form 

Ru = An = au ss + bu S Q + cuee + du s = 0 (4.20) 

wliere ci,b,c and d are functions of s and 0 obtained from chain rule. Seek 
approximate solutions of the form 

M * 

«*(»,«) = £ Y. “mt T m (s) (-1.21) 

m=0 |fc|<yV 

where T m is the tilth Chebyshev polynomial of the first kind, defined by 

r m (cos 0 ) = cos (m9) (4.22) 

and u N is determined by at the (M + 1) x 2 N Chebyshev-Fourier nodes 

(s t ,0j) = (cos iir/M,jir/N) (4.23) 

f = 0, M 

y = 0, 2N - 1 

At i = A/ the given Di rich let condition should be imposed. At i — 0,du/ds 
obtained from Ou/dr should be updated. Placing the inhomogeneous terms 
arising from the boundary conditions in a forcing term /, let us rewrite the 
final spectral operator as 

- f) {si,9j) = 0. . (4.24) 

The matrix formed by L, p a V is large and does not have sparse struc- 
tures. Thus in order to solve (4.24), an iterative procedure is desirable. We 
outline only the brief idea of implementing this procedure. Dropping the 
superscript N in (4.24) the method is as follows: 

^ n r 1 = + ^n(L sp U n f ) 
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(4.25) 



where a n is chosen so that the i 2 norm of the residual 


r n f l — / L 3p U n 


is minimum. This gives 

{ r m I J sp r n) 

( LspTn , pJ"n) 

where ( , ) denotes an f. 2 inner product. 


a n = - 


(4.26) 


(4.27) 


As a sample comparison we can generate a situation where the exact 
solution is 

u(x , </) = log ((x - -7) 2 + y 2 } (4.28) 

For the geometry of a circle (r = R(6) = 1) and for a .33 x 32 grid, C 2 
errors are listed for different values of R 0 0 . Solutions are compared against 
Bayliss, Gunzburger and Turkel’s procedure (BGT) with their first order 
(FO) and second order (SO) boundary conditions given in section two (Equa- 
tions 2.23). 10 CHE denotes the infinite order radiation condition of the 
present method. (Canuto, llariharan and Lustman). 


I t CO 

FO BGT 

SO BGT 

IO CHL 

1.2 

.16 x IQ -2 

.031 x 10" 2 

.00037 x 10 

3 

.15 x 10" 3 

.025 x 10~ 3 

.01 x 10" 5 

5 

.14 x 10 ~ 4 

.16 x 10~ 4 

.14 x 10" 4 


TABLE 1 

We see from this table that in all cases the infinite order radiation 
condition is superior, especially when the artificial boundary is near the body, 
which is desired anyway. As R 0 c increases the first order and second order 
conditions become better and comparable to the infinite order condition. 
This is because the grid size is fixed. If the grid size is increased, the infinite 
order condition will improve. Further illustrations can again be found in 
reference' 5 ! . 


Note further that this procedure can be easily extended to other types of 
problems. If we consider the Helmholtz equation A u+k 2 u = 0 with radiation 
condition at infinity we can write the outgoing solutions as follows: 

OO 

U (r,0)=£a„//< l >(fa-)e'"'' (4.30) 

n— 0 


where are the flankel functions of the first kind and order n. 


If we take the radial derivative of (4.30) we have 


du 


OO 


Q- = £ a n k H n l)> i kr ) e 


tnO 


(4.31) 


n~~ 0 


Again we note that a n tln\kr) is the nth Fourier coefficient of u(r, •), which 
can be inverted to give 

a„H^(kr) = i [ a(r, v )e-'"*dp (4.32) 

Substitution of this in (4.31) gives 

Ou ^ kHi i] '(kr) f 2 ’ , , 

= > — - / u(r,<p)e 

<> r Ill'Hkr) Jo 


p) d<p 


(1.33) 


or 


du 

= li o u 

dr 


(4.34) 


This is similar to what we had in equation (4.11) except that a finite number 
of values of llh^ and its derivatives must be calculated. To do this there are 
well known recurrence relations and numerical evaluations available in the 
literature. 


5. Boun dar y Con di tion s for T ime De pendent Problem s 

This has been the most difficult problem to handle numerically. It is 
still an open question if a nonlocal condition that is suitable for numerical 
calculations can be obtained. We saw in section 1 that it is relatively easy 
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to obtain an absorbing condition for a simple wave equation in one dimen- 
sion. We proceed in the same fashion to obtain boundary conditions in two 
dimensions. This process results in the procedure introduced by Engquist 
and Majdaf 7 l’l 8 L Before we proceed to discuss the two dimensional case, note 
that in three dimensions D’Alambert’s solution holds. If we consider 

^2 = u xx + ^yy + Uzz (5* 0 

then the general solution for spherically radiating and incoming waves can 
be written as 


u(r, t) — * f(r - ct) + - g(r + ct ) (5.2) 

r r 

This representation is well known and may be found in Morse and Feshback' 161 . 
The first part indicates outgoing waves. If one wants to impose a radiation 
condition, this part may be used to represent the solution as follows: 

u(r, t) = - f(r - ct) (5.3) 

r 

If we assume the sound speed is normalized to one, then it is easy to check 
from (5.3) that 

\l t + Ut -f- - =0 (5.-4) 

r 

which is local, both in space and in time, but provides an exact description 
of spherically radiating waves. 

Thus in three dimensions a radiation condition similar to what we dis- 
cussed in sections 2 through 5, in the time domain will become exact, if one 
imposes (5.1) at sulliciently far distances. In two as well as in three dimen- 
sions analogous to (5.2) one wants a cylindrically or spherically radiating 
waves at finite distances and that causes problems, as we will see a little 
later. 

Let us begin the discussion of two dimensional wave equations. Consider 
a wave traveling from left incident on the boundary x = L(L > 0) without 
reflecting (see figure 8) and governed by 




The operator B = d/dx + d/dt dictated the outgoing part of the wave 
given by ( 1.16). The operator d/dx — d/dt dictates the incoming part of the 
wave. To obtain a non reflective condition we set Bu = 0. Let us examine 
if it is possible to do an anologous argument for two dimensions. Equation 
(5.5) can be written as 



We have indicated by arrows as in (5.7) that in (5.8) the resulting opera- 
tors may have a similar meaning. Unfortunately we have a square root of an 
operator. But this can be easily explained using the theory of pseudo differ- 
entia! operators. For this purpose let us take Fourier transforms of equations 
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(5.5) with respect to t and y and call the corresponding dual variables r and 
This gives 

«xx + (r 2 - <; 2 )u = 0 (5.9) 

Equation (5.9) is a reduced wave equation in the x direction, whose solutions 
may be written as 

ape; r, f) = /t(r,c)e»N/ ra -^ x + B{t, ?) (5.10) 

Arrows in (5.10) indicate the waves moving in the right and left directions. 
In order to impose an outgoing solution we see that we must choose 


f r 2 ~ s * 2 x 


u ( x ; t ,$)= A(r,f)eV r 
Now inverting the Fourier transform in (5.11), we have 

«(*,»,<)= j I A(T,s)e-' r,+i '« H 'J ;r: ^*dTd< 


(5.11) 


(S- 12) 


Even though (5.12) gives a perfectly nonreflecting condition at r = L it 
is rather impractical to impose computationally. However, in the pseudo- 
diirerential operator terminology is/r 2 - <; 2 is nothing but the symbol of the 

pseudo differential operator J Engquist and Majda proceeded 

with this guideline that the approximations of this symbol i\J r 2 — £ 2 yield 
a family of approximate boundary conditions. In particular, approximations 
around £ = 0 (this has the physical meaning of near normal incidence) gives 
the boundary conditions that we are seeking. For example, the zeroth order 
approximation of i\/r 2 — c 2 is 


i\J T 2 - <; 2 za ir 

Now dilferentiate (5.12) with respect to x. 

u x (x,y,t)= I J is/r 2 -^ 2 ,l(r, <;) e trt *-'<» drd<; 

Substitute (5.13) into (5.1 1) and evaluate at x = L to obtain 

j j ir. ,l(r. f ) L drd<; 


(5.13) 


(5.11) 


(5.15) 
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The right hand side of (5.15) is nothing but the Fourier transform of — ut. 
Thus (5.15) can be rewritten as 


\ dx dt ) 


u 


= 0 


x=L 


(5.16) 


Which is the one dimensional condition (1.15). Obviously this will not be 
accurate enough for high resolution to handle two dimensional effects. I5n- 
gquist & Majda proceeded to obtain higher order approximations of the 
symbol i\J t 2 - <; 2 . The next order Taylor approximation is 

i\/ t' 2 - ( 2 ^ i (^t - (5.17) 

Substitution of (5.17) in (5.M) yields 


»,(i,M) = y I i(r-^) A( ( ,r)e-"‘ +i ^ +i ^ r ‘-^ L drd ! (5.18) 

Differentiating (5.18) with respect to t we see that 

Wxt — Hu 2 u yy-i ^ ^ ^ (5.19) 

In this way higher order boundary conditions can be generated. One warning 
should be given that the higher order Taylor approximants of the symbol 
i\J r 2 - <; 2 do not always yield stable boundary conditions. The proof of 
this statement is difficult and will be found in reference' 7 !. Instead these 
authors proposed Fade 1 approximants and found out that they are stable 
for all approximate boundary conditions. The second Taylor approximation 
coincides with the first Fade' approximation and from physical reasoning 
both the boundary conditions (5.16) and (5.14) are stable. 

There is another independent work due to Reynolds' 20 ', parallel to the 
work of Fngquist and Majda which derives the boundary condition (5.19) 
as a special case. Moreover, this work describes in detail on reducing edge 
reflections. Recently, Keys’ 2 ' 1 proposed a new method again by decomposing 
the wave equation into incoming and outgoing components to obtain a family 
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of boundary conditions. In this work, he derived Engquist and Majda’s as 
well as Reynolds’s conditions as special cases. 


Now suppose cylindrically radiating boundary conditions are imposed. 
Then it is necessary to change the equation (5.5) into cylindrical polar co- 
ordinates: 

n lt — «rr + “X u f. 90 H u r . (5.20) 

H r 

This has nonconstant coeIRcients and needs sufficient modification to apply 
the above theory. Again variable coefficient theory pseudo differential op- 
erators can be used. It is difficult to summarize this procedure, however 
there is a similar approach to that which we have given above. First Pade 
approximation of the resulting symbol yields the boundary condition 


( d d_ J_ \ 

\ dr dt 2 r ) 


u — 0 


to be imposed at sufficiently large distance r. 


( 5 . 21 ) 


It is interesting to note that this condition can be obtained from separa- 
tion of variables analogous to that of the spherically outgoing solution (5.3). 
In this case we do not have a simple form of D’Alembert’s solution, rather 
it has the form 


H(f , e , l)m mzll ( ao(e) + “jM + 

V r \ r 


(5.22) 


Eliminating <io(0) in (5.22) we see that 

(^ + I + i)' l = 0 ( 1/r5/2 )- 


(5.23) 


We see that boundary condition (5.21) agrees with the physical one (5.23). 
Eliminating U[(0),a 2 (0) etc., Bayliss & Turkel |2 l obtained a family of radi- 
ation conditions, which agree only with (5.23) of Engquist & Majda. Other 
higher order conditions differ from each other. Higher order conditions are 
very appealing, but difficult to implement. It will be of interest to both 
applied mathematicians as well as engineers to see how effective arc these 
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boundary conditions. Hariharan and Bayliss^ 10 ! implemented three dimen- 
sional version of (5.23) i.e., (5.1) to a practical problem described below. 
In three dimensions this boundary condition is asymptotically accurate to 

o(£). 

The problem is to solve for sound radiation into atmosphere from a 
cylindrical pipe. There is an incident wave on the left end of the pipe and 
sound radiates into atmosphere from open end of the pipe (figure 9). 


Inflow Boundary 


Far field Boundary 


Directivity Measurement 


Semi -Infinite Inlet 


\c. 




Figure 9 Computational plane of sound radiation problem 


The incident pressure wave, in particular has the form 

p(r,0,e)^f(r,e)e‘”‘* (5.24) 

where m is called the mode number. The governing wave equations have the 
form 


Op v 4- imw 

,, + u» + v r + — 0 

at r 

(5-25) 

'(- + ~ p =0 ■ 
Ot Oz 

(5.28) 

i)v dp 

at + Tr = 0 

(5.27) 

^ + ™ P = o 

at r 

(5.28) 


These equations are derived from Euler equations of the associated 
How problem and are available in reference 110 !. Due to cylindrical symmetry 
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(which one can recognize from these equations) it is sufficient to solve the 
problem in one plane. In this plane as in figure 5.2 the truncated boundaries 
are 71, 72 and 73 together with the axis of the pipe and the inflow boundary, 
where, a source term is to be imposed. 

There are boundary conditions on the inflow boundary, axis, and on 
the walls of the pipe. I3ut we shall be concerned only with the radiation 
condition on 71,72 and 73 . For rest of the details, the reader is referred to 
reference 1 10 l again. When the condition (5.23) is imposed for p, the acoustic 
pressure, we have 



?+^ = o 

dt 2 r 



Suppose a point on 7,(1 = 1,2,3) is at a distance R = \fz 2 + r 2 from the 
origin and the line from the origin to the point makes an angle a with the 
axis. Then 


dp 

OR 


dp dp , 

= cos a -(- — sin a 
dz dr 


Hut from (5.26) and (5.27) we see that 


(5..-S0) 


dp du 

dz dt 

dp dv 

dr dt 

Thus (5.20) becomes 


dp 

dt 


/ r 

(u cos a + v sin o) + - = 0 


(5.30) 


This condition, together with other boundary conditions, was used to 
solve the system (5.25) through (5.28) by a fourth order finite difference 
scheme to obtain solutions reported in figure 5.32. In this figure the vertical 
axis measures the sound pressure levels (dB) and the horizontal axis gives 
angles 0 where the calculations were made at a distance of 10 diameters 
of the pipe. For this situation when the inflow has a time dependence of 
the form e~ twt , Weiner- 1 lop f solutions can be computed. We compared our 
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procedure with the work of Savkar and Edelfelt^ 17 ! which uses the Weiner- 
Ilopf technique. Also, we compared the solutions with some experimental 
data from Viilc and Silcoxl t9 l The comparison is shown in figure 5.3 for a 
wave number k = 3.37 and for azimuthal angular dependence of the solution 
e l20 . Weiner- 1 lopf theory and the numerical solutions agree well. However, 
there is a discrepancy with the experimental results, especially near the axis 
for small values of 0. The reason is due to certain uncontrollable factors, 
such as plane wave and lower order mode of e t0 dependence, which arise in 
the experimental situations. The main point of emphasis here is that the 
radiation condition (5.23) is a suitable condition for a practical problem. 



0 10 20 30 40 50 60 70 80 90 

Angle, deg 


Figure 10 Comparison of results 
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